For purposes of illustration, the function Exp.cov.simple is
provided as a simple example and implements the R code discussed below.
It can also serve as a template for creating new covariance functions for the
Krig and mKrig functions. Also see the higher level function
stationary.cov to mix and match different covariance shapes and
distance functions. Functional Form: If x1 and x2 are matrices where nrow(x1)=m and
nrow(x2)=n then this function will return a mXn matrix where the (i,j)
element is the covariance between the locations x1[i,] and x2[j,]. The
covariance is found as exp( -(D.ij **p)) where D.ij is the Euclidean
distance between x1[i,] and x2[j,] but having first been scaled by theta.
Specifically if theta is a matrix to represent a linear
transformation of the coordinates, then let
u= x1%*% t(solve( theta)) and v= x2%*% t(solve(theta)).
Form the mXn distance matrix with elements:
D[i,j] = sqrt( sum( ( u[i,] - v[j,])**2 ) ).
and the cross covariance matrix is found by exp(-D).
The tapered form (ignoring scaling parameters) is a matrix with i,j entry
exp(-D[i,j])*T(D[i,j]). With T being a positive definite tapering function that is also
assumed to be zero beyond 1.
Note that if theta is a scalar then this defines an isotropic covariance
function and the functional form is essentially exp(-D/theta).
Implementation:
The function r.dist is a useful FIELDS function that finds
the cross Euclidean distance matrix (D defined above) for two sets of
locations. Thus in compact R code we have
exp(-rdist(u, v)**p)
Note that this function must also support two other kinds of calls:
If marginal is TRUE then just the diagonal elements are returned
(in R code diag( exp(-rdist(u,u)**p) )).
If C is passed then the returned value is
exp(-rdist(u, v)**p) %*% C
Radial basis functions Rad.cov: The functional form is
Constant* rdist(u, v)**p for odd dimensions
and Constant* rdist(u,v)**p * log( rdist(u,v)
For an m th order thin plate spline in d dimensions p= 2*m-d and must
be positive. The constant, depending on m and d, is coded in the fields
function radbas.constant. This form is only a generalized
covariance function -- it is only positive definite when restricted to
linear subspace. See Rad.simple.cov for a coding of the radial
basis functions in R code.
Stationary covariance stationary.cov:
Here the computation is to apply the function
Covariance to the distances found by the Distance function.
For example
Exp.cov(x1,x2, theta=MyTheta)
and
stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Exponential")
are the same. This also the same as
stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Matern",smoothness=.5).
Stationary tapered covariance stationary.taper.cov: The resulting
cross covariance is the direct or Shure product of the tapering function
and the covariance. In R code given location matrices, x1 and
x2 and using Euclidean distance.
Covariance(rdist( x1, x2))*Taper( rdist( x1, x2))
By convention, the Taper function is assumed to be identically
zero outside the interval [0,1]. Some efficiency is introduced within
the function to search for pairs of locations that are nonzero with
respect to the Taper. This search may find more nonzero pairs than
dimensioned by max.points. Given this error just pass a larger
for max.points explicitly. For spam.format TRUE the
multiplication with the C argument is done with the spam sparse
multiplication routines through the "overloading" of the %*%
operator. Currently this function only supports the Euclidean distance
function.
About the FORTRAN: The actual function Exp.cov and
Rad.cov calls FORTRAN to
make the evaluation more efficient this is especially important when the
C argument is supplied. So unfortunately the actual production code in
Exp.cov is not as crisp as the R code sketched above. See
Rad.simple.cov for a R coding of the radial basis functions.